# Updated table 1

# Function for manually adding stars
get_stars <- function(pval, ...){
  
  p <- pval
  
  star <- ifelse(
    p >= 0.1, "", 
    ifelse(
      p < 0.1 & p >= 0.05, "ASTRSK",
      ifelse(
        p < 0.05 & p >= 0.01, "ASTRSKASTRSK",
        ifelse(
          p < 0.01, "ASTRSKASTRSKASTRSK", NA)
      )
    )
  )
  
  return(star)
}


######################################################################################################################################################

# Prepare dataset
regression_vars <- dplyr::select(abba, treatment, uid, pweight, strata, block_code, dplyr::contains("y0"), fps_uid, uidcount)

bl_bal <- regression_vars %>%
  dplyr::left_join(., bl_dealer, by = "fps_uid") %>% 
  dplyr::left_join(
    ., el_dealer %>% 
      dplyr::mutate(block_code = as.character(block_code)) %>% 
      dplyr::select(-pweight, -block_code, -treatment), by = "fps_uid") %>% 
  dplyr::left_join(
    ., allocation_blk %>% 
      dplyr::select(treatment, strata, block_code, district_code, dis_value_perRC_total_y0,
                    dis_value_perRC_kero_y0, dis_value_perRC_rice_y0, dis_value_perRC_wheat_y0, 
                    dis_value_perRC_sugar_y0, dis_value_perRC_salt_y0)) %>% 
  dplyr::mutate(leakage = dis_value_perRC_total_y0 - value_total_y0,
                leakage_kero = dis_value_perRC_kero_y0 - value_kero_y0,
                leakage_rice = dis_value_perRC_rice_y0 - value_rice_y0,
                leakage_wheat = dis_value_perRC_wheat_y0 - value_wheat_y0,
                leakage_sugar = dis_value_perRC_sugar_y0 - value_sugar_y0,
                leakage_salt = dis_value_perRC_salt_y0 - value_salt_y0)

#####
lapply(c(1:length(c("leakage", "leakage_kero", "leakage_rice", "leakage_wheat", "leakage_sugar", "leakage_salt"))), function(x){ 
  form <- formula(paste0(c("leakage", "leakage_kero", "leakage_rice", "leakage_wheat", "leakage_sugar", "leakage_salt")[x], "~ treatment | strata | 0 | block_code"))
  reg  <- lfe::felm(form, weights = bl_bal$pweight, data = bl_bal) 
  reg$wald <- lfe::waldtest(reg, ~ treatment)[["p.F"]]
  return(c(c("leakage", "leakage_kero", "leakage_rice", "leakage_wheat", "leakage_sugar", "leakage_salt")[x], round(reg$beta, 2), round(reg$wald, 2)))
})

bl_bal %>% 
  dplyr::select(dplyr::all_of(c("leakage", "leakage_kero", "leakage_rice", "leakage_wheat", "leakage_sugar", "leakage_salt")), treatment, pweight) %>%
  dplyr::group_by(treatment) %>% 
  dplyr::summarise(dplyr::across(dplyr::all_of(c("leakage", "leakage_kero", "leakage_rice", "leakage_wheat", "leakage_sugar", "leakage_salt")), ~ weighted.mean(., w = pweight, na.rm = TRUE))) %>% 
  t(.) %>% 
  as.data.frame(.) %>% 
  janitor::row_to_names(row_number = 1) %>% 
  tibble::rownames_to_column(var = "Outcome") -> L

#####

# Create variables
bl_bal$b5_transaction_cost_y0 <- (bl_bal$b5_storing_ration_y0 + bl_bal$b5_transport_y0)/bl_bal$a14_rc_count
bl_bal$seeded_ind <- ifelse(bl_bal$uidcount > 0 & !is.na(bl_bal$uidcount), 1, NA)
bl_bal$seeded_ind <- ifelse(bl_bal$uidcount == 0 & !is.na(bl_bal$uidcount), 0, bl_bal$seeded_ind)
bl_bal %<>% dplyr::mutate(dplyr::across(dplyr::everything(), ~ ifelse(is.infinite(.), NA, .)))

# Check distributions
table(bl_bal$seeded_ind, bl_bal$treatment)
sum(is.na(bl_bal$seeded_ind[bl_bal$treatment == 1]))
sum(is.na(bl_bal$seeded_ind[bl_bal$treatment == 0]))

# Vector of variables for testing balance
variables_balance <- c(
  "value_total_y0", "value_rice_y0", "value_wheat_y0", "value_salt_y0", "value_sugar_y0", "value_kero_y0", 
  #"leakage",
  "c_cost_to_access_y0",
  #"seeded_ind",
  "b5_transaction_cost_y0"#,
  #"g2_caste_STSC_y0", "g2_caste_y0", "g6_land_sq_m_y0", "g7_ann_total_income_y0"
)

# Treatment and control means
means <- bl_bal %>% 
  dplyr::select(dplyr::all_of(variables_balance), treatment, pweight) %>%
  dplyr::group_by(treatment) %>% 
  dplyr::summarise(dplyr::across(dplyr::all_of(variables_balance), ~ weighted.mean(., w = pweight, na.rm = TRUE))) %>% 
  t(.) %>% 
  as.data.frame(.) %>% 
  janitor::row_to_names(row_number = 1) %>% 
  tibble::rownames_to_column(var = "Outcome")

# Regression adjusted difference and p-value
out    <- lapply(c(1:length(variables_balance)), function(x){ 
  form <- formula(paste0(variables_balance[x], "~ treatment | strata | 0 | block_code"))
  reg  <- lfe::felm(form, weights = bl_bal$pweight, data = bl_bal) 
  reg$wald <- lfe::waldtest(reg, ~ treatment)[["p.F"]]
  return(c(variables_balance[x], round(reg$beta, 2), round(reg$wald, 2)))
})

out_dat <- out %>% 
  as.data.frame(.) %>%
  t(.) %>% 
  tibble::as_tibble(.) %>%
  dplyr::rename(Outcome = V1, reg_diff = V2, pval = V3) %>% 
  dplyr::right_join(., means, by = "Outcome") 

out_dat <- cbind(out_dat, star = get_stars(pval = out_dat$pval)) %>% 
  dplyr::mutate(reg_diff = paste0(reg_diff, star)) %>% 
  dplyr::select(-star)

star_prep <- out_dat %>% 
  dplyr::select(Outcome, `1`, `0`, reg_diff, pval) %>% 
  dplyr::rename(Treatment = `1`, Control = `0`, `Regression-adjusted difference` = reg_diff) %>%
  dplyr::mutate(dplyr::across(dplyr::all_of(c("Treatment", "Control")), ~ round(as.numeric(.), 2)),
                Outcome = dplyr::case_when(
                  Outcome == "value_total_y0" ~ "Total value received",
                  Outcome == "value_rice_y0" ~ "Value received (rice)",
                  Outcome == "value_wheat_y0" ~ "Value received (wheat)",
                  Outcome == "value_salt_y0" ~ "Value received (salt)",
                  Outcome == "value_sugar_y0" ~ "Value received (sugar)",
                  Outcome == "value_kero_y0" ~ "Value received (kerosene)",
                  Outcome == "leakage" ~ "Leakage",
                  Outcome == "c_cost_to_access_y0" ~ "Beneficiary transaction costs",
                  Outcome == "seeded_ind" ~ "At least one Aadhaar seeded",
                  Outcome == "g2_caste_y0" ~ "Caste",
                  Outcome == "g2_caste_STSC_y0" ~ "SC/ST",
                  Outcome == "g6_land_sq_m_y0" ~ "Land owned (sq. meters)",
                  Outcome == "g7_ann_total_income_y0" ~ "Annual income",
                  Outcome == "b5_transaction_cost_y0" ~ "Dealer transaction costs",
                )) %>%
  as.matrix(.)

######################################################################################################################################################
###### Replicate original table 1 panels

### Blocks at baseline 
block_vars <- c(
  "ph_hh", "aay_hh", "uid_entry_hh", "offtake_hh_ph_rice", "offtake_hh_aay_rice", "fps_block", "median_fam_size"
  #, "secc_percent", "newdata_percent", "nodata_percent"
)

block_means <- block_rand_char %>% 
  dplyr::select(dplyr::all_of(block_vars), treatment) %>%
  dplyr::group_by(treatment) %>% 
  dplyr::summarise(dplyr::across(dplyr::all_of(block_vars), ~ mean(., na.rm = TRUE))) %>% 
  t(.) %>% 
  as.data.frame(.) %>% 
  janitor::row_to_names(row_number = 1) %>% 
  tibble::rownames_to_column(var = "Outcome")

block_out    <- lapply(c(1:length(block_vars)), function(x){ 
  form <- formula(paste0(block_vars[x], "~ treatment | strata | 0 | block_code"))
  reg  <- lfe::felm(form, data = block_rand_char) 
  reg$wald <- lfe::waldtest(reg, ~ treatment)[["p.F"]]
  return(c(block_vars[x], round(reg$beta, 2), round(reg$wald, 2)))
})

block_out_dat <- block_out %>% 
  as.data.frame(.) %>%
  t(.) %>% 
  tibble::as_tibble(.) %>%
  dplyr::rename(Outcome = V1, reg_diff = V2, pval = V3) %>% 
  dplyr::right_join(., block_means, by = "Outcome") 

block_out_dat <- cbind(block_out_dat, star = get_stars(pval = block_out_dat$pval)) %>% 
  dplyr::mutate(reg_diff = paste0(reg_diff, star)) %>% 
  dplyr::select(-star)

block_star_prep <- block_out_dat %>% 
  dplyr::select(Outcome, `1`, `0`, reg_diff, pval) %>% 
  dplyr::rename(Treatment = `1`, Control = `0`, `Regression-adjusted difference` = reg_diff) %>%
  dplyr::mutate(dplyr::across(dplyr::all_of(c("Treatment", "Control")), ~ round(as.numeric(.), 2)),
                Outcome = dplyr::case_when(
                  Outcome == "ph_hh" ~ "Priority HH",
                  Outcome == "aay_hh" ~ "AAY HH",
                  Outcome == "uid_entry_hh" ~ "Aadhaard numbers seeded per RC",
                  Outcome == "offtake_hh_ph_rice" ~ "Rice disbursed per priority HH",
                  Outcome == "offtake_hh_aay_rice" ~ "Rice disbursed per AAY HH",
                  Outcome == "fps_block" ~ "Number of FPS",
                  Outcome == "median_fam_size" ~ "Median HH size",
                  Outcome == "secc_percent" ~ "% of rationcard holders identified via SECC",
                  Outcome == "newdata_percent" ~ "% of rationcard holders identified by application",
                  Outcome == "nodata_percent" ~ "% of rationcard holders without eligibility info")) %>%
  as.matrix(.)

### Aadhaar numbers at baseline 
aadhaar_vars <- c("seeded_ind", "seeded_missing")

aadhaar_means <- aadhaarstatus %>% 
  dplyr::select(dplyr::all_of(aadhaar_vars), treatment) %>%
  dplyr::filter(!is.na(treatment)) %>%
  dplyr::group_by(treatment) %>% 
  dplyr::summarise(dplyr::across(dplyr::all_of(aadhaar_vars), ~ mean(., na.rm = TRUE))) %>% 
  t(.) %>% 
  as.data.frame(.) %>% 
  janitor::row_to_names(row_number = 1) %>% 
  tibble::rownames_to_column(var = "Outcome")

aadhaar_out    <- lapply(c(1:length(aadhaar_vars)), function(x){ 
  form <- formula(paste0(aadhaar_vars[x], "~ treatment | strata | 0 | block_code"))
  reg  <- lfe::felm(form, data = aadhaarstatus) 
  reg$wald <- lfe::waldtest(reg, ~ treatment)[["p.F"]]
  return(c(aadhaar_vars[x], round(reg$beta, 2), round(reg$wald, 2)))
})

aadhaar_out_dat <- aadhaar_out %>% 
  as.data.frame(.) %>%
  t(.) %>% 
  tibble::as_tibble(.) %>%
  dplyr::rename(Outcome = V1, reg_diff = V2, pval = V3) %>% 
  dplyr::right_join(., aadhaar_means, by = "Outcome") 

aadhaar_out_dat <- cbind(aadhaar_out_dat, star = get_stars(pval = aadhaar_out_dat$pval)) %>% 
  dplyr::mutate(reg_diff = paste0(reg_diff, star)) %>% 
  dplyr::select(-star)

aadhaar_star_prep <- aadhaar_out_dat %>% 
  dplyr::select(Outcome, `1`, `0`, reg_diff, pval) %>% 
  dplyr::rename(Treatment = `1`, Control = `0`, `Regression-adjusted difference` = reg_diff) %>%
  dplyr::mutate(dplyr::across(dplyr::all_of(c("Treatment", "Control")), ~ round(as.numeric(.), 2)),
                Outcome = dplyr::case_when(
                  Outcome == "seeded_ind" ~ "% of HH with at least one Aadhaar seeded",
                  Outcome == "seeded_missing" ~ "% of HH missing any Aadhaar seeded")) %>%
  as.matrix(.)

### First stage
fs_vars <- c("e1_dealer_has_epos", "e2_dealer_used_epos_3", "e2_dealer_used_epos_2", "e2_dealer_used_epos_1")

fs_means <- abba %>% 
  dplyr::select(dplyr::all_of(fs_vars), treatment, pweight) %>%
  dplyr::filter(!is.na(treatment)) %>%
  dplyr::group_by(treatment) %>% 
  dplyr::summarise(dplyr::across(dplyr::all_of(fs_vars), ~ weighted.mean(., w = pweight, na.rm = TRUE))) %>% 
  t(.) %>% 
  as.data.frame(.) %>% 
  janitor::row_to_names(row_number = 1) %>% 
  tibble::rownames_to_column(var = "Outcome")

fs_out    <- lapply(c(1:length(fs_vars)), function(x){ 
  form <- formula(paste0(fs_vars[x], "~ treatment | strata | 0 | block_code"))
  reg  <- lfe::felm(form, data = abba) 
  reg$wald <- lfe::waldtest(reg, ~ treatment)[["p.F"]]
  return(c(fs_vars[x], round(reg$beta, 2), round(reg$wald, 2)))
})

fs_out_dat <- fs_out %>% 
  as.data.frame(.) %>%
  t(.) %>% 
  tibble::as_tibble(.) %>%
  dplyr::rename(Outcome = V1, reg_diff = V2, pval = V3) %>% 
  dplyr::right_join(., fs_means, by = "Outcome") 

fs_out_dat <- cbind(fs_out_dat, star = get_stars(pval = fs_out_dat$pval)) %>% 
  dplyr::mutate(reg_diff = paste0(reg_diff, star)) %>% 
  dplyr::select(-star)

fs_star_prep <- fs_out_dat %>% 
  dplyr::select(Outcome, `1`, `0`, reg_diff, pval) %>% 
  dplyr::rename(Treatment = `1`, Control = `0`, `Regression-adjusted difference` = reg_diff) %>%
  dplyr::mutate(dplyr::across(dplyr::all_of(c("Treatment", "Control")), ~ round(as.numeric(.), 2)),
                Outcome = dplyr::case_when(
                  Outcome == "e1_dealer_has_epos" ~ "Dealer has an ePOS machine at endline survey",
                  Outcome == "e2_dealer_used_epos_3" ~ "Dealer used ePOS in January 2017",
                  Outcome == "e2_dealer_used_epos_2" ~ "Dealer used ePOS in February 2017",
                  Outcome == "e2_dealer_used_epos_1" ~ "Dealer used ePOS in March 2017")) %>%
  as.matrix(.)

bal_tab_original <- rbind(block_star_prep, aadhaar_star_prep, star_prep, fs_star_prep)

balance_table_final <- 
  stargazer::stargazer(
    bal_tab_original,
    type = "latex", 
    align = FALSE, 
    digits = 2,
    column.sep.width = "5pt",
    df = FALSE, 
    digits.extra = 2,
    nobs = TRUE,
    float = FALSE,
    header = FALSE,
    model.numbers = TRUE,
    multicolumn = FALSE,
    dep.var.caption = "",
    title = "", 
    label = "tab:balance",
    notes.append = FALSE, 
    notes.align = "l",
    notes = ""
  )

# Make title better 
balance_table_final[5] <- "\\multicolumn{1}{l}{} & \\multicolumn{1}{c}{Treatment} & \\multicolumn{1}{c}{Control} & 
                            \\multicolumn{1}{c}{\\specialCellCenter{Regression-\\\\adjusted \\\\ difference}} & \\multicolumn{1}{c}{$p$-value}\\\\"  
balance_table_final <- 
  starpolishr::star_insert_row(
    balance_table_final, 
    c(" \\cmidrule(lr){2-2}  \\cmidrule(lr){3-3}  \\cmidrule(lr){4-4}  \\cmidrule(lr){5-5}"),
    insert.after = c(5)
  )

balance_table_final <- 
  starpolishr::star_insert_row(
    balance_table_final, 
    c("& \\multicolumn{1}{c}{(1)} & \\multicolumn{1}{c}{(2)} & \\multicolumn{1}{c}{(3)} & \\multicolumn{1}{c}{(4)}\\\\"),
    insert.after = c(6)
  )

# Align first column to the left
balance_table_final[2] <- "\\begin{tabular}{@{\\extracolsep{5pt}} lcccc}"

# Remove leading zeros
balance_table_final <- stringr::str_replace_all(string = balance_table_final, pattern = "0\\.", replace = ".")

balance_table_final

# Make asterisks better
balance_table_final <- gsub(x = balance_table_final, pattern = "ASTRSK", replacement = "*")

# Add panel separators
balance_table_final <- 
  starpolishr::star_insert_row(
    balance_table_final, 
    c("\\multicolumn{5}{@{} l}{\\textit{ Panel A: Baseline characteristics (n = $132$ blocks)}} \\\\ \\addlinespace"),
    insert.after = c(8)
  )

balance_table_final <- 
  starpolishr::star_insert_row(
    balance_table_final, 
    c("\\bottomrule \\addlinespace \\multicolumn{5}{@{} l}{\\textit{ Panel B: Survey balance (n = $3960$ households)}} \\\\ \\addlinespace"),
    insert.after = c(18)
  )

balance_table_final <- 
  starpolishr::star_insert_row(
    balance_table_final, 
    c("\\bottomrule \\addlinespace \\multicolumn{5}{@{} l}{\\textit{ Panel C: Program implementation (n = $3578$ ration cards)}} \\\\ \\addlinespace"),
    insert.after = c(27)
  )

# Save
starpolishr::star_tex_write(
  starlist = balance_table_final,
  file = paste0(here::here(), "/replication/exhibits/TableX_updated_survey_balance_table.tex"), 
  headers = FALSE
)

######################################################################################################################################################
###### investigate leakage

allocation_blk %>% 
  dplyr::select(block_code, district_code, treatment, strata, dplyr::contains("dis_value_perRC") & !dplyr::contains("imp")) %>% 
  tidyr::pivot_longer(-c(block_code, district_code, treatment, strata), names_pattern = "(.*)_(.*)", names_to = c(".value", "item")) %>%
  dplyr::group_by(treatment, item) %>%
  dplyr::summarise(dplyr::across(dplyr::contains("dis_value"), ~ mean(., na.rm = TRUE))) %>% 
  tidyr::pivot_longer(cols = dplyr::contains("dis_value")) -> A

abba %>% 
  dplyr::select(block_code, district_code, treatment, strata, dplyr::contains("value_") & !dplyr::matches("ent_|stat_ent_|_paid_|_fps_")) %>% 
  tidyr::pivot_longer(-c(block_code, district_code, treatment, strata), names_pattern = "(.*)_(.*)", names_to = c(".value", "item")) %>%
  dplyr::group_by(treatment, item) %>%
  dplyr::summarise(dplyr::across(dplyr::contains("value_"), ~ mean(., na.rm = TRUE))) %>% 
  dplyr::select(!dplyr::matches("mar|apr|may")) %>%
  tidyr::pivot_longer(cols = dplyr::contains("value_")) -> B 

C <- rbind(A, B)

D <- cbind(A, B) %>% 
  dplyr::mutate(value = `value...4` - `value...8`) %>% 
  dplyr::select(treatment = `treatment...1`, item = `item...6`,  name = `name...7`, value) %>% 
  dplyr::mutate(name = paste0(name, "_leakage"))

E <- rbind(C, D)

ggplot() + 
  geom_line(
    dplyr::filter(
      E, stringr::str_detect(string = name, pattern = "dis_") == TRUE),
      mapping = aes(
        x = factor(item, levels = c("y0", "jan17", "feb17", "mar17")), 
        y = value, color = name, group = name)) + 
  geom_point(E, mapping = aes(x = factor(item, levels = c("y0", "jan17", "feb17", "mar17")), y = value, color = name, group = name)) + 
  scale_color_manual(values = c(
    #"dis_value_perRC_kero" = "grey50",
    "dis_value_perRC_rice" = "blue",
    #"dis_value_perRC_salt" = "grey50",
    #"dis_value_perRC_sugar" = "grey50",
    #"dis_value_perRC_wheat" = "grey50",
    "dis_value_perRC_total" = "red",
    #"value_kero" = "grey50",
    "value_rice" = "steelblue",
    #"value_salt" = "grey50",
    #"value_sugar" = "grey50",
    #"value_wheat" = "grey50",
    "value_total" = "lightcoral",
    #"value_kero_leakage" = "grey50",
    "value_rice_leakage" = "forestgreen",
    #"value_salt_leakage" = "grey50",
    #"value_sugar_leakage" = "grey50",
    #"value_wheat_leakage" = "grey50",
    "value_total_leakage" = "palegreen"),
    labels = c(
      #"dis_value_perRC_kero" = "Kerosene (disbursed)",
      "dis_value_perRC_rice" = "Rice (disbursed)",
      #"dis_value_perRC_salt" = "Salt (disbursed)",
      #"dis_value_perRC_sugar" = "Sugar (disbursed)",
      #"dis_value_perRC_wheat" = "Wheat (disbursed)",
      "dis_value_perRC_total" = "Total (disbursed)",
      #"value_kero" = "Kerosene (received)",
      "value_rice" = "Rice (received)",
      #"value_salt" = "Salt (received)",
      #"value_sugar" = "Sugar (received)",
      #"value_wheat" = "Wheat (received)",
      "value_total" = "Total (received)",
      #"value_kero_leakage" = "Kerosene (leakage)",
      "value_rice_leakage" = "Rice (leakage)",
      #"value_salt_leakage" = "Salt (leakage)",
      #"value_sugar_leakage" = "Sugar (leakage)",
      #"value_wheat_leakage" = "Wheat (leakage)",
      "value_total_leakage" = "Total (leakage)")) + 
  labs(x = "Month", y = "Value", color = "") +
  facet_wrap(. ~ treatment, labeller = as_labeller(c("0" = "Control", "1" = "Treatment"))) + 
  theme_bw() + 
  theme(text = element_text(size = 14))

# Can't reject that t and c values are different for value disbursed or received. 
t.test(allocation_blk$dis_value_perRC_rice_y0[allocation_blk$treatment == 1], allocation_blk$dis_value_perRC_rice_y0[allocation_blk$treatment == 0])
t.test(abba$value_rice_y0[abba$treatment == 1], abba$value_rice_y0[abba$treatment == 0])

lfe::felm(dis_value_perRC_rice_y0 ~ treatment | strata | 0 | block_code, data = allocation_blk) %>% summary(.)
lfe::felm(value_rice_y0 ~ treatment | strata | 0 | block_code, weights = abba$pweight, data = abba) %>% summary(.)

t.test(bl_bal$leakage[bl_bal$treatment == 1], bl_bal$leakage[bl_bal$treatment == 0])
t.test(bl_bal$leakage_kero[bl_bal$treatment == 1], bl_bal$leakage_kero[bl_bal$treatment == 0])
t.test(bl_bal$leakage_rice[bl_bal$treatment == 1], bl_bal$leakage_rice[bl_bal$treatment == 0])
t.test(bl_bal$leakage_wheat[bl_bal$treatment == 1], bl_bal$leakage_wheat[bl_bal$treatment == 0])
t.test(bl_bal$leakage_sugar[bl_bal$treatment == 1], bl_bal$leakage_sugar[bl_bal$treatment == 0])
t.test(bl_bal$leakage_salt[bl_bal$treatment == 1], bl_bal$leakage_salt[bl_bal$treatment == 0])

lfe::felm(leakage ~ treatment | strata | 0 | block_code, data = bl_bal, weights = bl_bal$pweight) %>% summary(.)
lfe::felm(leakage_kero ~ treatment | strata | 0 | block_code, data = bl_bal, weights = bl_bal$pweight) %>% summary(.)
lfe::felm(leakage_rice ~ treatment | strata | 0 | block_code, data = bl_bal, weights = bl_bal$pweight) %>% summary(.)
lfe::felm(leakage_wheat ~ treatment | strata | 0 | block_code, data = bl_bal, weights = bl_bal$pweight) %>% summary(.)
lfe::felm(leakage_sugar ~ treatment | strata | 0 | block_code, data = bl_bal, weights = bl_bal$pweight) %>% summary(.)
lfe::felm(leakage_salt ~ treatment | strata | 0 | block_code, data = bl_bal, weights = bl_bal$pweight) %>% summary(.)

